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Abstract We consider integrated circuits with semiconductors modeled by modi- 
fied nodal analysis and drift-diffusion equations. The drift-diffusion equations are 
discretized in space using mixed finite element method. This discretization yields a 
high dimensional differential-algebraic equation. We show how proper orthogonal 
decomposition (POD) can be used to reduce the dimension of the model. We com- 
pare reduced and fine models and give numerical results for a basic network with 
one diode. Furthermore we discuss an adaptive approach to construct POD models 
which are valid over certain parameter ranges. Finally, numerical investigations for 
the reduction of a 4-diode rectifier network are presented, which clearly indicate 
that POD model reduction delivers surrogate models for the diodes involved, which 
depend on the position of the semiconductor in the network. 
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1 Introduction 

In this article we investigate a POD-based model order reduction for semiconduc- 
tors in electrical networks. Electrical networks can be efficiently modeled by a 
differential-algebraic equation (DAE) which is obtained from modified nodal anal- 
ysis. Denoting by e the node potentials and by and jy the currents of inductive 
and voltage source branches, the DAE reads (see ||8l [T9l ) 
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-A//.v(f), (1) 

0, (2) 
vs(t). (3) 

Here, the incidence matrix A — [AR,Ac,Ai,Ay ,Ai] represents the network topology, 
e.g. at each non mass node /, Oij = 1 if the branch j leaves node / and = — 1 if 
the branch j enters node /. If node / and branch j are not connected, then = 0. 
qc, g and (^l are continuously differentiable functions defining the voltage-current 
relations of the network components. The continuous functions and 4 are the 
voltage and current sources. For a basic example consider the network in Figure[T] 
Under the assumption that the Jacobians 

Dc{e,t):=^{e,t), DG{e,t) := ^{e,t), DlUj) ^U,t) 
ae ae a j 

are positive definite, analytical properties (e.g. the index) of DAE ([B-© are inves- 
tigated in fTTl and [ 14 1. In linear networks, the matrices Dq, Dq and Dl are positive 
definite diagonal matrices with capacitances, conductivities and inductances on the 
diagonal. 



Ac—qc{Aj:e,t) +A«^(Aje,f) +ALiL+Aviv = 

j^<t>L{iL.t)-Ale^ 
Ale = 



Fig. 1 Basic test circuit witli 
one diode. Tlie network is 
described by 



Av = ( 1, 
A« = ( 0, 



0)\ 



ei{t) hit) 



sit) 







n{t,x) 
p{t,x) 



e2{t) 



Often semiconductors themselves are modeled by electrical networks. These 
models are stored in a library and are stamped into the surrounding network in order 
to create a complete model of the integrated circuit. Here we use a different approach 
which models semiconductors by the transient drift-diffusion equations. Advantages 
of this approach are the higher accuracy of the model and fewer model parameters. 
On the other hand, numerical simulations are more expensive. For a comprehensive 
overview of the drift-diffusion equations we refer to ||Tl|2]|4l|T0l[T5l. Using the nota- 
tion introduced there, we have the following system of partial differential equations 
for the electrostatic potential \//(f , jc), the electron and hole concentrations n{t,x) and 
p{t,x), and the current densities Jn{t,x) and Jp{t,x): 
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qdtn +div/„ 
qd,p + div Jp 



q{n-p-C), 

qR{n,p,J,„Jp), 
-qR(n,p,Jn,Jp), 

^pqi-Ur'^p-pVxt/), 



with (f ,x) e [0, T] X £2 and n CR'' (d = 1 ,2,3). The nonhnear function R describes 
the rate of electron/hole recombination, q is the elementary charge, e the dielectric- 
ity, /x„ and jXp are the mobihties of electrons and holes. The temperature is assumed 
to be constant which leads to a constant thermal voltage Ut- The function C is the 
time independent doping profile. 

This paper is organized as follows. In Section |2] we present the (scaled) mathe- 
matical model for the complete coupled network including semiconductors modeled 
by the DD equations. In Section [3] we describe the numerical treatment of the net- 
work using the method of lines. For the spacial discretization of the DD models we 
use mixed finite elements. For the time integration of the resulting DAE system we 
use DASSL with step width control and so obtain snapshots y''{ti, •), i — 1, . . . ,k, 
which represent the state of the circuit and the semiconductors at time . Based on 
these snapshots and POD we construct a reduced model in Section ID We discuss 
the properties of the reduced model with respect to parameter changes in Section|5] 
In order to obtain a reduced model which is valid over a considerable range of pa- 
rameters it is necessary to refine the model. For this purpose we adapt the reduced 
basis method combined with the greedy approach proposed in [iZJ to our setting. 
In Section|6]we finally present numerical experiments, and also discuss advantages 
and shortcomings of our approach. 



2 Complete coupled system 

In the present section we develop the complete coupled system for a network with 
ns semiconductors. We will not specify an extra index for semiconductors, but we 
keep in mind that all semiconductor equations and coupling conditions need to be 
introduced for each semiconductor. 

For the sake of simplicity we assume that to a semiconductor m semiconductor 
interfaces Fg^k ^ F C d Q , k = I , . . . ,m ctre associated, which are all Ohmic contacts, 
compare Figure |2] The dielectricity e shall be constant over the whole domain Q . 
We focus on the Shockley-Read-Hall recombination 



R{n,p): 



np — 7]^ 



Zp{n + ri) + z„{p + ri) 
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which does not depend on the current densities. Herein, t„ and Xp are the average 
hfetimes of electrons and holes and rj is the constant intrinsic concenti'ation which 
satisfy Tj^ —np if the semiconductor is in thermal equilibrium. 

The coupling of the drift-diffusion equations and the electrical network is estab- 
lished as follows. The current A^y'^ is added to equation ([T]i, where 



js,k^ / {J„+Jp-ed,V\j/)-vda. (4) 

I.e. the current is the integral over the current density J„ + Jp plus the displacement 
current in normal direction v. Furthermore, the node potentials of nodes which are 
connected to a semiconductor interface are introduced in the boundary conditions 
of the drift-diffusion equations (see also Figure|2]i: 



(5) 



«(f,x) = -( VC(x)2+4tj2 + C(x)), (6) 



x)2 +4772 -C(x) , (7) 



for {t,x) E [0,T] X Fo^k- Here, \lfhi{x) denotes the build-in potential and rj the con- 
stant intrinsic concentration. All other parts of the boundary are isolation boundaries 
I] := r \ To, where Vy/ • V = 0, 7„ • V = and 7p • v = holds. 

The complete model forms a partial differential-algebraic equation (PDAE). The 
analytical and numerical analysis of such systems is subject to current research, see 
in |7] [17] [l9l . The simulation of the complete coupled system is expensive and nu- 
merically difficult due to bad scaling of the drift-diffusion equations. The numerical 
issues can be significantly reduced by the unit scaling procedure discussed in ifTSl . 
That means we substitute 

x = Lx, xi/ = UtY, «=||C|1„«, p = \\C\\^p, C=|1C|1„C, 

Jfl - jlfjJfj^ Jp - jXpJp^ TJ )7||C||cxa, 



where L denotes a specific length of the semiconductor The scaled drift-diffusion 
equations then read 
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XA\l/ = n-p-C, (8) 

-<9,n + V„divy„ = Rin,p), (9) 

d,p + Vpdi\Jp = -R{n,p), (10) 

7„= Vn-nV\j/, (11) 

Jp-^-yp-p'^W, (12) 



where we omit the tilde for the scaled variables. The constants are given by A := 
ITMaZ^ v„ and Vp 



Fig. 2 Sketch of a coupled 
system with one semiconduc 
tor. Here 

for all {t,x) e [0,r] xFo,!. 



3 Simulation of the full system 

Classical approaches for the simulation of drift-diffusion equations (e.g. Gummel 
iterations |6 1) approximate J„ and Jp by piecewise constant functions and then solve 
equations (|9]l and ( fTOl i with respect to n and p explicitly. This helps reducing the 
computational effort and increases the numerical stability. For the model order re- 
duction approach proposed in the present work this method has the disadvantage of 
introducing additional non-linearities, arising from the exponential structure of the 
Slotboom variables, see [17|. 

Since the electrical field represented by the (negative) gradient of the electrical 
potential plays a dominant role in (l8Tl-(fT2]i and is present also in the coupling con- 
dition dU, we provide for it the additional variable g^f — Vy/ leading to the following 
mixed formulation of the DD equations: 
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Adivgi^ = n-/5-C, (13) 

-<9,n + V„divy„ = R{n,p), (14) 

d,p + VpdiwJp = -R{n,p), (15) 

8v= Vy/, (16) 

J„= yn-ngxf,, (17) 

Jp^~'^P~Pgv (18) 

The weak formulation of ([T3ll-([T8l) then reads: Find \i/,n,p e [0, T] x L^{Q) and 
gif,,J„Jp e [0,T] X //o.Af(div,i2) such that 



A / div^^(p= / {n-p)(p- C(p, (19) 

Ji2 Ji2 Ji2 

d,n(p + V„ div7„ ^ = / R{n,p) (p, (20) 



i2 

•JrP (P + Vp div/p (p = - R{n,p) (p, (21) 

i2 Ji^ Ji2 



/ 



8v ■(!> = - V/div0+/ (22) 

J„-(j) = - n div<j)+ n<j)-V- ng«,-<j), (23) 
12 Jn Jr Jn 



Jp-(p= pdiw(p- p(p-v- pg^-(p, (24) 
n Jn Jr Jn 

are satisfied for all (p £ L?{£2) and (j) e //o^A'(div,i2) where the space HQi^{div,£2) 
is defined by 

i/(div,i2) {y e L^(r2)'' : divy eL^{Q)}, 
//o,iv(div,i2) {y e //(div,f2) : v^Oonf?}. 

Consequently, the boundary integrals on the right hand sides in equations (|22]|-(|24]| 
reduce to integrals over the interfaces Fp.k, where the values of y/, n and p are 
determined by the Dirichlet boundary conditions (|5]l-(|7]l. We note that, in contrast to 
the standard weak form associated with ©-(O, the Dirichlet boundary values are 
naturally included in the weak formulation (fT9]l-(l24b and the Neumann boundary 
conditions have to be included in the space definitions. This is advantageous in 
the context of POD model order reduction since the non-homogeneous boundary 
conditions (|5]l-(|7]i are not present in the space definitions. 

Here, equations (fT9]l- (l24l i are discretized in space with Raviart-Thomas finite 
elements of degree (RTq), alternative discretization schemes for the mixed problem 
are presented in |4|. To describe the TJTo-approach for d = 2 spatial dimensions, let 

be a triangulation of £2 and let S" be the set of all edges. Let (oj :— {E £ S' : E C 
/}} be the set of edges at the isolation (Neumann) boundaries. The potential and the 
concentrations are approximated in space by piecewise constant functions 
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/(f),«''(f),/(f) e Lh {y e L^n) : y\T{x) = cj, W e 
and the discrete fluxes g'^{t), J'!,{t) and Jp{t) are elements of the space 

RTq:= {y: Q. ^W' : y\T{x) ^ + bTX, 07- e M'', /^r e R, 

blfi • V£ =0, for all inner edges £}. 

Here, \y]E denotes the jump 7^ —yW- over a shared edge E of the elements and 
T . The continuity assumption yields RTq C //(div,i2). 
We set 

i//,,o,iv(div,i2) := (RTon//o,;v(div,f2)) C //o.7v(div,i2). 

Then it can be shown, that ///,.o.7v posses an edge-oriented basis {^j} j=\....M- We 
use the following finite element Ansatz in (fT9]l-(l24]l 



N 



M 



(25) 



i=l .7=1 

A' M 

n''(f,x) = £n,-(f)(PiW, -/n^f = E JnMiix), 

1=1 ,7=1 

W M 

'■=1 i=i 

where := |^|, i.e. the number of elements of and M := l^?! — {S'nI, i.e. the 
number of inner and Dirichlet boundary edges. 
This in ([T9]l-(|24li yields 

M r- N C C 



N 



n 



' E "'(0 / (Pi (Pk + V„y J„j{t) / div (j)j (pk - / R{n'',p'') (pk = 0, 
~i Jn Jn Jn 



i=l 



a 



a 



a 



M 



N 



L / 0r ^/ + E V'Kf ) / 9/ div = w'' <l>rv 



M 



7=1 
M 



E / 0r ^/ + E "'(0 / div 0/ + / nV^, • 0/ = n''<l)rV 



E ■^Pj(0 / 0; • 0/ - E/'-^O / <P, div(/)/ + / •(/), = -//(/),• V, 



i2 
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which represents a nonUnear, large and sparse DAE for the approximation of the 
functions n, p, gxj/, J„, and Jp. In matrix notation it reads 



/ \ / 

-Mz.n(f) 

MLp{t) 









+ 



V / 



-Ml Ml XD 



VpD 



M„ 



Mh 



\ 



MhJ 



n{t) 
Pit) 

Jn{t) 



with 



Afem 



( \ 

- S^R{n\p^^) Cf, 
JnR{n\p'')(p 




b:-- 



+ :^{n\p\^^)=b{e{t)), 



( -InCcp \ 



V -Irp'<P-y J 



and 



/ 

Jn 



R{n\p'')<P :-- 



( SnR{n\p'')<pi\ 



All other integrals in and h are defined analogously. The matrices Ml € M^^^ 

and Mh G m.^^^ are mass matrices in the spaces L/, and Hf,fi^N^ respectively, and 
D g M'VxM jj^g gjj^ jjQ^ jgjjg^ jjjg fpjjj^ 



Residual Based Sampling in POD Model Order Reduction 

Problem 1 (full model). 

Acj^qc{Ale{t),t) +A«g(Aje(f ),f ) +Az,iL(f) +Ay;V(0 



/ \ 

-Mifiit) 

MLp{t) 







^FEM 



+Asjs{t)+Aih{t)=Q, (26) 

j^ML{t).t)-Ale{t)=Q, (27) 

ATe(r)-v,(0=0, (28) 
js{t)-C,Jn{t)-C2Jp{t)-Cig^{t)=Q, (29) 

n(t) 

g% +^(«''/'4)-Me(0)=0, (30) 

Jn{t) 
\jp{t)l 



where ( |29] l represents the discretized linear coupling condition (|4]l. 



4 Model reduction 

We now aim to reduce the computational effort of repeated dynamical simulations 
by applying proper orthogonal decomposition (POD) to the drift-diffusion equa- 
tions. The idea is to replace the large number of local model-independent ansatz 
and test functions {0,}, {^y} by only a few nonlocal model-dependent Ansatz func- 
tions for the respective variables. 

The snapshot variant of POD introduced in fTSl works as follows. We run a sim- 
ulation of the unreduced system and collect / snapshots ij/'^itk, •), n^'{tk, •), p^{tk, 
8%{h,-)^ j'nifk,-), Jpih,-) at time instances f^. e {fi,...,f/} C [0,r]. The optimal 
selection of the time instances is not considered here. We use the time instances 
delivered by the DAE integrator 

Since every component of the state vector y := {y/,n,p,gxi/,Jn,Jp) has its own 
physical meaning we apply POD MOR to each component separately. Among other 
things this approach has the advantage of yielding a block-dense model and the 
approximation quality of each component is adapted individually. 

Let X denote a Hilbert space and let y'' : [0, T] x X — > W with some r e N. 
The Galerkin formulation ( |25] l yields y'^{t,-) S Xi, := span{0^, . . . }, where 
{(^f}i<j<n denote n linearly independent elements of X. The idea of POD consists 
in finding a basis {m' , . . . , m"'} of the span of the snapshots 

span 1 /(f,, •) = f^yf'ipfi-), with ^ = 1, . . . ,/ I 
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satisfying 

{u\...,u'}= argmin £ /{tk,-) -t,{y\h,-)y{-))xv\- ' 

{v',...,v'}cX k=\ 1=1 

for 1 < s < m, where I <m<l. The functions {u'}i<i<s are orthonormal in X and 
can be obtained with the help of SVD as follows. 

Let the matrix Y := {y'''^ ,■ ■ ■ € W^' contain as columns the coefficient vec- 
tors of the snapshots. Furthermore, let M := {{^f ■^J)x)i<ij<n be the positive def- 
inite mass matrix with its Cholesky factorization M = LLJ . Let {U ,L,V) denote 
the singular value decomposition of Y := L^Y, i.e. Y = UZV'^ with U G W", 
V e M'^', and a matrix E € R"^' containing the singular values Ci > CJ2 > . . . > 
Cm > Oin+i > ■ • • > CT/ > 0. We Set U := L"^f7(. i-^y Then, the s-dimensional POD 
basis is given by 

span|M'(-) = J = l,...,i|. 

The information content of {m',. . . ,m*} with respect to the scalar product {■,-)x with 



0<A{s) = ^^^^<l, (31) 

is given by 1 — A{s). Here 4(5') measures the lack of information of {m\. ..,«•'} 
with respect to span{>'''(fi ,■),... ,y^{ti,-)}. 

The POD basis functions are now used as trial and test functions in the Galerkin 
method. 

If the snapshots satisfy inhomogeneous Dirichlet boundary conditions POD is 
performed for 

f{t) = \l/{t)-\l/r{t), h{t) = n{t)-nr{t), p{t)=p{t)-pr{t), 

with (//, , Tir, Pr denoting reference functions satisfying the Dirichlet boundary con- 
ditions required for xj/, n and p. This guarantees that the POD basis admits homoge- 
neous boundary conditions on the Dirichlet boundary. In the case of the mixed finite 
element approach the introduction of a reference state is not necessary, since the 
boundary values are included more naturally through the variational formulation. 

The time-snapshot POD procedure delivers Galerkin Ansatz spaces for n, p, 
gY, Jn and Jp. This leads to the Ansatz 



(32) 
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The injection matrices 



nNxsu 



Uj„ e R^'''-'- 



Up el 



dMxsj 



contain the (time independent) POD basis functions, the vectors 7(.) the correspond- 
ing time-variant coefficients. The numbers are the respective number of POD 
basis functions included. AssembUng the POD system yields the DAE 



/ \ 

-%{t) 

Ut) 




V ) 



iPOD 



with 



7n{t) 
YJnit) 



ApoD — U AfEuU 

I 



V 



-U^MlU„ U^MlU,, XU^DU,^ 



\ 



I 



and U = diag{Uti/,U„,Up,Ug^,Uj„,Ujp). Note that we exploit the orthogonaHty of 
the POD basis functions, e.g. uJMlU„ = l/jMiUp = Inxn and Uj^MnUg^ = 

UJMhUji, ~ UJMhUjp = Imxm- The arguments of the nonlinear functional have 
to be interpreted as functions in space. 

All matrix-matrix multiplications are calculated in an offline-phase. The nonlin- 
ear functional ^ has to be evaluated online. A possible reduction method for the 
nonlinearity is called Discrete Empirical Interpolation and is discussed in |5j. Here, 
we focus on the speed-up in solving linear systems in the implicit time integrator, 
which are now small. The reduced model for the network now reads 
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Problem 2 (reduced order model). 



Acj^qc{Ale{t),t) +Affg(Aje(f ),f ) +Az,jz,(f) +Ay;V(0 

+Asjs{t)+A,U{t)^0, (33) 

j^ML{t),t)-Ale{t)^Q, (34) 

ATe(f)-y,(f) = 0, (35) 

js (f ) - C, f/y„ 7y„ (f ) - C2f/y„ 7y„ (0 - CiU,, (0 = 0, (36) 



-U\nn'''''y'"',^^'"')-U^b{e(t))=Q. (37) 



/ \ 












7p(f) 



+ ApoD 


Yp{t) 
7«„(f) 







7j„(0 


V j 







5 Residual-based sampling 



We now consider parameter dependent models. Possible parameters are physical 
constants of the semiconductors (e.g. length, permeability, doping) and parameters 
of the network elements (e.g. frequency of sinusoidal voltage sources, value of re- 
sistances). We do not distinguish between inputs and parameters of the model. 

Let there be r € N parameters and let the space of considered parameters be given 
as a bounded set ^ C W. We construct the reduced model based on snapshots from 
a simulation at a reference parameter (Oi G One expects that the reduced model 
approximates the unreduced model well in a small neighborhood of coi , but one 
cannot expect that the reduced model is valid over the complete parameter set 
In order to create a suitable reduced order model we consider additional snapshots 
which are obtained from simulations at parameters 02,03,... G The iterative 
selection of (0^+1 at a step k is called parameter sampling. Let P/^ denote the set of 
selected reference parameters, P<. := {oi , O2, . . . , (Oi^} C 

We neglect the discretization error of the finite element method and its influence 
on the coupled network and define the error of the reduced model as 

^(o;P):=z''(o)-z™^(o;P), (38) 

where z''(o) :— (e''(o), 7^(0), 7^(0), 3'''(o))^ is the solution of problem[T]at the 
parameter o with discretized semiconductor variables 37'' : = ( 1/'' , n'' , p'' , g'^ , jj^ iJpY ■ 
z^^^{(i);P) denotes the solution of the coupled system in problem |2] with reduced 
semiconductors, where the reduced model is created based on simulations at the 
reference parameters P C The error is considered in the space X with norm 
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\\z\\x:^ {WehMjvhMkh, 

II VllL2([0,r],L2(i2))' Il"llz.2([0,r],i2(i2)) JpllL2([0,r]X2(i2)): 
ll^vllL2{[0,7'],//o,A,(div,i2))' 

ll-^«llL2([0,r].//o.Jv(div.i2))' ll-^pllL2([0,7-],//o.jv(div,i2))^ 

Obvious extensions apply when there is more than one semiconductor present 
Furthermore we define the residual ^ by evaluation of the unreduced model 
(l30b at the solution of the reduced model z^^"{(0;P), i.e. 



POD 



\ 







'S:(„POD POD POD\ 



^FEM 



)-b{e''0O(t)), (39) 



Note that the residual of equations (I26]l-(|29]l vanishes. 

We note that the same definitions are used in |9| for linear descriptor systems. 
In 1(9 1 an error estimate is obtained by deriving a linear ODE for the error and ex- 
ploiting explicit solution formulas. Here we have a nonlinear DAE and at the present 
state we are not able to provide an upper bound for the error ||ff (a);P)||x which 
would yield a rigorous sampling method using for example the Greedy algorithm 
of d. 

We propose to consider the residual as an estimate for the error The evaluation 
of the residual is cheap since it only requires the solution of the reduced system and 
its evaluation in the unreduced DAE. It is therefore possible to evaluate the residual 
at a large set of test parameters Ptest C Similar to the Greedy algorithm of |[T2| . 
we add to the set of reference parameters the parameter where the residual becomes 
maximal. 

The magnitude of the components in error and residual may be large and a proper 
scaling should be applied. For the error we consider the component-wise relative 
error, i.e. 

111/(60) -v/^Q^(a);P)||^2([o,r],^2(^)) ||«^'(co)-n^o^(co;P)||^2([o.r],^2(^)) 

llr''(«)llL2([0,r],L2(i2)) ' l|n''(«)llL2([o,r],L2(i2)) 

and the residual is scaled by a block-diagonal matrix containing the weights 
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\ 



d„{(0)I 



dp{(0)I 



\ 



dj„{(0)lj 



The weights c/(.)(a)) > may be parameter-dependent. These weights are chosen 
in a way that the norm of the residual and the relative error are component-wise 
equal at the reference frequencies (Ok where we know 7!^ {(Ok) from simulation of the 
unreduced model, i.e. 



ing the weights df^.){(0\), d(^.\^{(Ok) piecewise linearly. Extrapolation is done by 
nearest-neighbour interpolation to ensure the positivity of the weights. 
We summarize our ideas in the following sampling algorithm: 

Algorithm 1 (Sampling) 

1. Select (0\ e Ptest C 3^, tol > 0, and set k 1, Pi := {(0\}. 

2. Simulate the unreduced model at Oi and calculate the reduced model with POD 
basis functions U\. 

3. Calculate weight functions t/(.)(a)) > according to i40i for all (Ok G Pk- 

4. Calculate the scaled residual ||D(co).^(z^'^^(co,7\))|| for all (0 G Ptest- 

5. Check termination conditions, e.g. 



• no progress in weighted residual. 

6. Calculate (Ok+\ := arg max^^^^^^^ ||D((o)^(z™^((0,f<.))||. 

7. Simulate the unreduced model at (Ok+i and create a new reduced model with 
POD basis Uk+i using also the already available information at (0\, (Ok. 

8. Set Pk+\ U {o^^+i}, k:=k-\-\ and goto 3. 

The step 7 in Algorithm [T] can be executed in different ways. If offline time and 
offline memory requirements are not critical one may combine snapshots from all 
simulations of the full model and redo the model order reduction on the large snap- 
shot ensemble. Otherwise we can create a new reduced model at reference frequency 
(Ok+i with POD-basis U and then performing an additional POD step on {Uk,U)- 




H</(a),)-v/^^^(a),;P)H^2(p,rix2(i2)) 




. maxa,eP,„, ||Z)(«)^(z™«(tO,P,))|l < tol, 
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The discussed finite element method is implemented for spatial one-dimensional 
(ID) problems in MATLAB. The resulting high dimensional DAE is integrated us- 
ing the DASSL software package 1 13 1. We assume that the differentiation index of 
the network is 1. Otherwise one should switch to alternative integrators. In order to 
solve the Newton systems efficiently which arise from the BDF method, the vari- 
ables of the sparse system are reordered with respect to minimal bandwidth. 

A basic test circuit with a single diode is depicted in Figure [1] where the model 
parameters are presented in Table [1] The input Vs{t) is chosen to be sinusoidal with 
amplitude 5 [V]. The numerical results in Figure[3]show the capacitive effect of the 
diode for high input frequencies. Similar results are obtained in [16] using the sim- 
ulator MFCS. In the sequel the frequency will be considered as a model parameter. 



Table 1 Diode model parameters. 



Parameter 


Value 




Parameter 


Value 




L (length) 


lo-'* 


[cm] 


e 


1.03545- 10-'2 


[F 1 cm] 


Ut 


0.0259 


[V] 


ri 


1.4-10'° 


Wlcm"] 




1350 


[cm^ / (V sec)] 




330-10-" 


\sec] 




480 


[cm^/lvsec^] 




33-10-" 


\sec] 


a (contact area) 


10-5 


[cm^] 


C{x), x<L/2 


-9.94 -10'5 


[l/cm^] 








C(x), x>L/2 


4.06- 10'* 


Wlcm'] 


frequency 


= 1 MHz 


frequency = 1 GHz 


frequency = 


5 GHz 




0.5 1 0.5 1 0.5 1 

scaled time scaled time scaied time 



Fig. 3 Current through the basic network for input frequencies 1 MHz, 1 GHz and 5 GHz. The 
capacitive effect is clearly demonstrated. 
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Fig. 4 Relative error between 
reduced and unreduced prob- 
lem at the fixed frequency 
IQio [Hz\. 




Fig. 5 Time consumption for 
simulation runs of Figure |4] 
The dashed line indicates 
the time consumption for the 
simulation of the original full 
system. 



18 
16 

„ 14 

o 

i£. 12 


.i 10 



- reduced simulation 
full simulation 




A(S) 



We first validate tlie reduced model at a fixed reference frequency of 10^" [Hz\. 
Figure |4| shows the development of the relative error between the reduced and the 
unreduced numerical solutions, plotted over the lack of information Zi , see dSTT i. The 
number of POD basis functions for each variable is chosen such that the indicated 
approximation quality is reached, i.e. A := 4,^ ~ 4,, ~ Zip ~ A^,^ — ^J,, — ^ • Since 
we compute all POD basis functions anyway, this procedure does not involve any 
additional costs. 

In Figure |5]the simulation times are plotted versus the lack of information A. 

Figure |6] shows the total number of singular vectors s required to guarantee a 
certain information content \ — A. It can be seen that with the number of singu- 
lar vectors included increasing only linearly, the lack of information tends to zero 
exponentially. 

We now apply Algorithm [T] to provide a reduced order model of the basic cir- 
cuit and we choose the frequency of the input voltage as model parameter As 
parameter space we chose the interval := [10^, lO'^] [Hz\- We start the investi- 
gation with a reduced model which is created from the simulation of the full model 
at the reference frequency Oi := lO'" [Hz\ ■ The number of POD basis functions s is 
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450 



Fig. 6 The number of re- 
quired singular values grows 
only logarithmically with 
the requested information 
content. 




chosen such that the lack of information A{s) is approximately 10~'. The relative 
error and the weighted residual are plotted in Figure|7] We observe that the weighted 
residual is a rough estimate for the relative approximation error. Using Algorithm[T] 
the next additional reference frequency is a>2 10^ [Hz] since it maximizes the 
weighted residual. The second reduced model is constructed on the same lack of 
information A :~ 10^^. Here we note that in practical applications, the error is not 
known over the whole parameter space. 



sampling step 1 



Fig. 7 Relative reduction er- 
ror (solid line) and weighted 
residual (dashed line) plotted 
over the frequency parameter 
space. The reduced model is 
created based on simulations 
at the reference frequency 
ffli := 10'" [H7]. The refer- 
ence frequencies are marked 
by vertical dotted lines. 



10" 



10" 



10" 



10 



10 



- error 
residual 
reference frequencies 




10 

parameter (frequency) 



10 



The next two iterations of the sampling algorithm are depicted in Figures[8]and|9l 
Based on the residual in step 2, one selects 0)3 := 1.0608 • 10^ [Hz] as the next 
reference frequency. Since no further progress of the weighted residual is achieved 
in step 3, the algorithm terminates. The maximal errors and residuals are given in 
Table |2] 
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Fig. 8 Relative reduction er- 
ror (solid line) and weighted 
residual (dashed line) plotted 
over the frequency parameter 
space. The reduced model is 
created based on simulations 
at the reference frequen- 
cies coi := IQi" [Hz\ and 
e)2 := 10** [//z]. The reference 
frequencies are marked by 
vertical dotted lines. 



10 



10 



10" 



10 



10 



lO" 



10' 
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sampling step 2 



- error 
residual 
reference frequencies 




10 



10'" 10'' 

parameter (frequency) 



Fig. 9 Relative reduction er- 
ror (solid line) and weighted 
residual (dashed line) plotted 
over the frequency parameter 
space. The reduced model 
is created based on simu- 
lations at the reference fre- 
quency coi := lO'" [Hz\, 
a>2 := 10** [Hz], and 
e>3 := 1.0608-10'' [Hz]. 
The reference frequencies 
are marked by vertical dotted 
lines. 



10' 
10' 
10° 
10' 
10"' 
10"^ 
10- 



sampling step 3 



- error 
residual 
reference frequencies 




10" 



10'" 10'' 

parameter (frequency) 



Table 2 Progress of refinement method. 

step k reference parameters max. scaled residual max. relative error 

(at frequency) (at frequency) 

1 {1.0000- 10'°} 9.9864-102 3.2189-10" 

(1.0000-10**) (1.0000-10**) 



{1.0000-10*, 


1.5982-10-2 


4.3567-10-2 


1.0000-10'°} 


(1.0608-10') 


(3.4551-10') 


{1.0000- 10^ 


2.2829 - 10-2 


1.6225-10-2 


1.0608- 10^ 


(2.7283-10') 


(1.8047-10"") 


1.0000 - lO'O} 
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Finally we note that the presented reduction method accounts for the position 
of the semiconductors in a given network in that it provides reduced order models 
which for identical semiconductors may be different depending on the location of 
the semiconductors in the network. The POD basis functions of two identical semi- 
conductors may be different due to their different operating states. To demonstrate 
this fact, we consider the rectifier network in Figure [TO] Simulation results are plot- 
ted in Figure [TTl The distance between the spaces and f/2 which are spanned, 
e.g., by the POD-functions t/^ of the diode and t/^ of the diode ^2 respectively, 
is measured by 

d{U^,U^):— max min ||m — v||2. 

[|«I|2=1 lk||2 = l 

Exploiting the orthonormality of the bases f/^ and and using a Lagrange frame- 
work, we find 

d{u\u^) = \J2-2Vx, 

where X is the smallest eigenvalue of the positive definite matrix SS^ with Sij = 
(m^ j)2. The distances for the rectifier network are given in Table[3] While the 
reduced model for the diodes Si and ^3 are almost equal, the models for the diodes 
Si and ^2 are significantly different. Similar results are obtained for the reduction of 
n, p, etc. 



eiit) 



Fig. 10 Rectifier network. 




Table 3 Distances between reduced models in the rectifier network. 



A 


d{U\U-) 




10-* 


0.61288 


5.373- 10-* 




0.50766 


4.712- 10-* 


10 * 


0.45492 


2.767-10-'' 


lo-'? 


0.54834 


1.211-10-*' 
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Fig. 11 Simulation results 
for the rectifier network. The 
input Vj is sinusoidal with 
frequency 1 [GHz] and offset 
+1.5 [V]. 
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frequency = 1 GHz 
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